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Abstract 

Utilizing the Haar transform, we study the higher order spectral properties of mean field 
avalanche models, whose avalanche dynamics are described by Poisson statistics at a critical point 
or critical depinning transition. The Haar transform allows us to obtain a time series of noise 
powers, H(fi,t), that gives improved time resolution over the Fourier transform. Using H(f\,t) 
we analytically calculate the Haar power spectrum, the real 1.5 spectra, the second spectra, and 
the real cross second spectra in mean field avalanche models. We verify our theoretical results 
with the numerical results from a simulation of the T = mean field nonequilibrium random field 
Ising model (RFIM). We also extend our higher order spectra calculation to data obtained from a 
numerical simulation of the T = infinite range RFIM for d = 3, and experimental data obtained 
from an amorphous alloy, Fe<i\Co§^B\^. We compare the results and obtain novel exponents. 
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I. INTRODUCTION 



There are disordered systems that respond to slow driving with discrete jumps or 
avalanches with a broad range of sizes, referred to as crackling Such crackling sys- 
tems are characterized by many interacting degrees of freedom and strong interactions that 
make thermal effects negligible, such as: charge density waves, vortices in type II supercon- 
ductors, crack propagation, earthquakes, and Barkhausen noise in magnets. The avalanche 
dynamics of the systems mentioned above are described by Poisson statistics (given by Eq. 
(JIJ) in mean field theory at a critical point or critical depinning transition [3, S|- 

Through our analysis of mean field avalanche models we determine the following spectral 
functions: Haar power spectrum, real 1.5 spectra, second spectra, and real cross second 
spectra for systems that have avalanche dynamics given by Eq. (0). This analysis pro- 
vides new tools for noise analysis in dynamical systems, since there are very few theoretical 
calculations of higher order noise statistics. 

The Haar transform allows us to obtain a power versus time series, H(f\,t), needed to 
calculate higher order spectra. These higher order spectra give valuable information about 
;he avalanche dynamics in Barkhausen noise not accessible through ordinary power spectra 
. Higher order spectra also have been used to obtain crucial information about a variety 



of diverse systems such as: metastable states in vortex flow p , natural auditory signals 7] , 
conductance-noise in amorphous semiconductors fluctuating current paths in devices 9|, 
and quasi-equilibrium dynamics of spin glasses While much experimental work has 



glasses iu|. wmie i 



been done studying higher order spectra |4|, |5j, |6|, [TJ, |8|, |9|, [10|, we present a rigorous mean 
field treatment that is applicable to a broad range of systems Q]. This analysis will allow 
a better understanding of the dynamics of these systems, and provide a direct method of 
comparison to experiment or observation. 

In addition, we compare our general results from mean field theory to Barkhausen noise 
obtained: from a mean field simulation of the T = random field Ising model (RFIM) ^j] , 
from a simulation of the T = infinite range RFIM (IRM) in d = 3, and from experiment. 
We also find novel exponents from our analysis, and we compare our results from theory, 
simulation, and experiment; we find key similarities and differences. 
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II. THE MODEL 



A. Mean Field RFIM 



The T = mean field RFIM consists of an array of iV spins (sj = ±1), which may 

point up (sj = +1) or down (sj = —1). Spins are coupled to all other spins (through 

a ferromagnetic exchange interaction J), and to an external field H(t) which is increased 

adiabatically slowly. To model disorder in the material, we assign a random field, hi, to 

each spin, chosen from a distribution P(hi) = exp(—hi/2R 2 )/\^2irR, where R determines 

the width of the Gaussian probability distribution and therefore gives a measure of the 

amount of quenched disorder for the system. The Hamiltonian for the system at a time t 

is given by: H = — X)( JM + H(t) + h^Si, where M = -k s j is the magnetization of the 

« i 
system. Initially, H(—oo) = — oo and all the spins are pointing down. Each spin is always 

aligned with its local effective field h\ = JM + H(t) + 



B. RFIM with infinite range Forces (IRM) in 3D 

The T = Infinite Range RFIM (IRM) consists of a 3D lattice of N spins, with a 
Hamiltonian at a time t given by: H = J2 —Js i Sj — J2i(H(t) + h i — J in fM)si, where J in f > 

<ij> 

is the strength of the infinite range demagnetizing field, and (ij) stands for nearest neighbor 
pairs of spins. The local effective field is given by h\ = JJ2<ij> s j + H(t) + hi — J in fM. 
The addition of a weak Ji n f ~ 4? to the traditional RFIM causes the system to exhibit 

rf n n 

self-organized criticality (SOC) |12t Il3l Il4j. This means that as H is increased the model 



always operates at the critical depinning point, and no parameters need to be tuned to 
exhibit critical scaling behavior (except ^ — > 0). We limit our analysis to a window of 



H values where the slope of M(H) is constant and the system disj 
behavior. Details of the simulation algorithm are given elsewhere [11 



ays front propagation 



C. Avalanche Dynamics in the T = RFIM 

The external field H(t) is adiabatically slowly increased from — oo until the local field, 
h e /* , of any spin Sj changes sign, causing the spin to flip P, [r|. It takes some microscopic 
time At for a spin to flip. The spin flip changes the local field of the coupled spins and may 
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cause them to flip as well, etc. This avalanche process continues until no more spin flips are 
triggered. Each step of the avalanche, that is each At, in which a set of spins simultaneously 
flip, is called a shell. The number of spins that flip in a shell is directly proportional to the 
voltage V(t) during the interval At that an experimentalist would measure in a pick-up coil 
wound around the sample. In our simulations we denote the number of spins flipped in a 
shell at a time t by n t (= V(t)). The first shell of an avalanche (one spin flip) is triggered by 
the external field H(t), while each subsequent shell within the avalanche is triggered only 
by the previous shell, since H(t) is kept constant while the avalanche is propagating. H(t) 
is only increased when the current avalanche has stopped, and is increased only until the 
next avalanche is triggered (i.e. ^ — > 0). The number of shells in an avalanche times At 
defines the pulse duration, T, or the time it took for the entire avalanche to flip. The time 
series of n t values for many successive avalanches creates a Barkhausen train analogous to 
experiment. 



III. EXPERIMENT 



We compare our results for theory and simulation with results obtained from an exper- 
iment performed on an (unstressed) amorphous alloy, Fe2\Co^B\^. Measurements were 
performed on a 21 cm x 1 cm x 30 fim ribbon of F e<i\C o^B^ alloy, a soft amorphous fer- 
romagnet obtained from Gianfranco Durin. The domain walls run parallel to the long axis 
of the material, with about 50 domains across the width. A solenoid, driven with a triangle 
wave, applies a magnetic field along the long axis of the sample. Since domain wall motion 
dominates over other means of magnetization in the linear region of the loop, data were col- 
lected in only a selected range of applied fields near the center of the loop. The Barkhausen 
noise was measured by a small pick-up coil wound around the center of the sample. This 
voltage signal was amplified, anti-alias filtered and digitized, with care taken to avoid pick- 
up from ambient fields. Barkhausen noise was collected for both increasing and decreasing 
fields for 80 cycles of the applied field through a saturation hysteresis loop. The driving 
frequency was 0.01 Hz; this corresponds to c = 0.09, where c is a dimensionless parameter 
proportional to the applied field rate and is defined in the Alessandro Beatrice Bertotti 
Montorsi model (ABBM model) for the Barkhausen effect in metals In this way, our 
measurements should be well inside the c < 1 regime identified in the ABBM model, in which 
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we can expect to find more or less separable avalanches rather than continuous domain wall 
motion. 



IV. THEORY 

A. Poisson Distribution 

The probability distribution for the avalanche dynamics in the class of mean field 
avalanche models we are interested in is given by 3, [J: 

P(n = l,n 1 ,n 2 ,...,n oo ) = r TT (1) 

em! fj~ 2 n t \ 

The above probability distribution (Eq. is for the time series of a single infinite 
avalanche at the critical point. In the context of the mean field RFIM, n t represents the 
number of spins flipped at a time t. That is, each n t represents a shell of the avalanche, and 
since an avalanche begins with one spin flip we have that n = 1. 

Let (.) represent the average over Eq. ^ In order to calculate the Haar transform and 
higher order spectra in mean field theory we need the following quantities, where m > I > 
k > and j > 0: 



(n,> = 1 (2) 

{n)) =3 + 1 (3) 

(njn j+k ) = (nf ) =j + l (4) 

K 3 ) = ^(3j 3 + 5j + 2) (5) 

K) = ^(6j 3 + 13j 2 + 9j + 2) (6) 

(n jnj+k n j+l ) = (nj) + k(n]) (7) 

( n j n j+k n j+l n j+m) 

= (up + (2k + /)<nj> + ~(k + 21- 1)(„J) (8) 



These relations can be determined from Eq. (^Q). We notice that the time ordering of 
the indices plays on important role. Details of how Eqs. (J2HHJ) were derived are given in 
Appendix |XJ 
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B. Haar Power 



The Haar transform is a simple wavelet transform with basis states consisting of single- 
cycle square waves. We use the Haar transform rather than the Fourier transform since it 
gives us improved time resolution in exchange for less frequency resolution. Time resolution 
is important to our purpose since we are interested in studying how the power contribution 
around a frequency f\ changes along the duration of the avalanche. 

Physically, the Haar power, H(t,fi), is the absolute square of the time integral over a 
period (of duration l//i) of a single-cycle square wave times a section of the train centered 
around t |4|. In other words, to determine the Haar power we integrate a square wave of 
duration l/fi over a section of the noise train centered at time t; this integrated segment 
is then squared to assure our resulting values is positive definite, this squared segment 
corresponds to H(t,fi), the Haar power at time t and around frequency f\. In order to 
analytically determine the Haar power series we first define the sum over N shells: 



N/2-1 

xt{N)= Y, nm+iAt (9) 
1=0 

N-l 

= Y n m+At (10) 

l=N/2 

Where t = Ni, and At is the time separation between shells. The Haar power series for 
a frequency f\ = is defined as: 

H(fi = J^ t ^ = ^) = i(*t(N) - x-(N)r) (11) 
To evaluate Eq. E]we determine the following relations: 



(xf(N)) = NAt/2 (12) 
N/2-1 

((xt(N)) 2 )= E [2(iV/2-l-m) + l]« + J (13) 

m=0 
N-l 

((xr(A0) 2 > = £ [2(N-l-m) + l)(n Ni+m ) (14) 

m=N/2 

N N/2 - 1 

(xf(N)x.(N)) = - £ (n Ni+m ) (15) 

/ m=0 
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When evaluating Eq. (jl3H15j) we must take into account the time ordering of the indices, 
since the time ordering is needed to evaluate the ensemble average of 2-pt (Eq. Hj), 3-pt (Eq. 
EJ), and 4-pt (Eq. |SJ) correlation functions. Refer to Appendix [B] for details. 

The above sums can be evaluated with the help of Eq. (J3j). Now using Eqs. (|12j) - (jl5j) we 
obtain the following exact result: 

m = i^,i) = k^-+h as) 



To find the Haar power spectrum we sum over all i, that corresponds to the sum over the 
Haar wavelets. To do this we define the maximum duration of the avalanche to be T. Now 
since blocks of N shells have been summed over, we perform a sum over Eq. (jlfijl from up 
to T/N - 1: 



i T/N-l 1 

S B (Jl) = 7r E #(/i = T7,0 (17) 

1 i=0 iv 

= ^(2(At) 2 + -^) (18) 
- — 2 (19) 



12/f 

The Haar power spectrum, Snifi), is in excellent agreement with the Haar power spec- 
trum determined from simulation. The Fourier power spectrum, Sp(fi) = \M differs 
by an additive constant and a constant factor from Sn{fi)- The additive constant, which 
depends on At, is in effect a result of aliasing produced by the discrete sampling. 

C. 1.5 Spectra, Second Spectra, and Cross Second Spectra 

The 1.5 spectra, second spectra, and cross second spectra are defined below: 



S, 5 (h, h) - m jjy t (20) 



s ^ h) ~ witjm (21) 



S2(/2 ' fb > fa) ~ (H(t a j a )) ta m b j b )) tb (22) 
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Where v(t, fx) is the sum of n« around the time t over a duration l/fi, and f 2 is the 
frequency conjugate to t. Also, F t {...} is the discrete Fourier transform with respect to 
t. Now (...) t denotes sum over t over the duration of the avalanche (J2t=o )> similarly 
(•••)*„ = £t=o and (-)ti = Et=o ■ For tli e cross second spectra, S 2 (f 2 , f b , f a ), we 
require that /& > / a . In addition, in the definition of S 2 (f 2 , fb, fa) we require two times, t a 
and tb. The time t a labels the starting point of a single-cycle square wave with a period 
of 1/ f a along the Barkhausen train, and therefore t a takes on values that are multiples of 
1/ f a - Similarly, t b labels the starting point of a single-cycle square wave with a period of 
1/ fb, and takes on values that are multiples of 1/ f b . 

To calculate Eq. (|2(JH22j) we first write the product of Fourier transforms as the Fourier 
transform of a convolution. This mathematical identity allows us to separate the Fourier 
transform (F t {■■■}) from the ensemble average ((.)) This leaves us with the following quan- 
titles: (v(t, f x )H{t + 6, (H(t, h)H{t + 6, f\)), and (H(t a + 6, f a )H(t b , f b )) where 9 is 
the convolution variable. These quantities may then be rewritten as a sum of 3-pt or 4-pt 
correlation functions, and subsequently evaluated. We obtain the general scaling forms: 



ri. 5 (/i) = T ~Y!\v(t, h)H(t, A)) ~ 4sn ( 23 ) 



t=0 Jl 

2 



A 

t=o Ji 



r 2 (A)= Yl (^«/i)>^i (24) 



T-l/fata + l/fa-l/fb JJ 3 

r 2 (/ 6 ,r)=r £ £ (H(t b J b )H(ta,f a ))~-^- (25) 

t a =0 tb—ta Jb 

Re{S 1 . 5 (f2ji)} = ^ + r[ N 5 \h) (26) 

J2 



B 



5 2 (/ 2 ,/i) = ^ + rf ) (/ 1 ) (27) 



Re{S 2 (f 2 , f a , fb)} = % + 4 N \f b , r) (28) 

J2 

The exponents V1.5, V 2 , Q1.5, Q2 are given in Table H]for mean field theory, the IRM, and 
the experiment. Also, Ai 5 , A 2 , -B1.5, and B 2 are non-universal constants. In Eq. (|25j) and 
Eq. (|2HD r = fa/h- Plots of Eqs. are given in Figs. ©-(ED- 

The functions Fi. 5 (/i), r 2 (/i), and T 2 (f b , r) are independent of f 2 since they resulted from 



the 6 = evaluation of (u(t, h)H(t+0, f x )), (H(t, fi)H{t+6, f{j), and (H{t a +6, f a )H(t b , f£ 
(see Fig. [TJ). These / 2 independent terms are referred to as Gaussian background terms 
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) 



Now ri. 5 ; (/i), and T ( 2 >(f b ,r) are r 1 . 5 (/ 1 ), r 2 (/i), and T 2 (f b ,r) normalized by 

(H(t,fx)) t , (H(t,fx))l, and (H(t b , f b )) tb (H(t a , f a )) ta , respectively. The first terms of Eq. 
(J2HJ) and Eq. (|27jl have no fi dependence, since the fi dependence of these terms drops 
out after they are normalized by (H(t,fi)) t and (H (t, fx))1 , respectively. Nevertheless, the 
lack of fi dependence in the first terms of second spectra and real 1.5 spectra is in excellent 
agreement with our simulation results, see Fig. 121 and the inset of Fig. El Also, we defer 
the discussion of Im{S' 2 (/2, fa, /&)} an d -^{Si.st^, /i)}> since they are very sensitive to the 
non-stationary properties of the infinite avalanche in the analytical model. Please refer to 
Appendix O for details of how Eqs. (j23H28j) are calculated. 

We present our results in Eq. (J23H28)) in terms of general scaling forms since the mean 
field and IRM simulations, as well as the experimental data, obey the same scaling form as 
mean field theory, only with different exponents and non-universal constants. 

V. SIMULATION 

A. Mean Field Simulation 

We perform 300 runs of a simulation of the mean field RFIM. We collect data taken from 
H e [0, 0.00125] at R = 0.79788 (R c = 0.79788456) in systems with iV = 15 x 10 6 spins, 
and J = 1 . From this data we determine our simulation results agree with the scaling forms 
given in Eq. (|23H28j) with exponents given by Table HI See Fig. EEl 

B. Infinite Range Model Simulation 

We perform 60 runs of a 3D simulation of the IRM. The data was taken from H 6 
[1.25, 1.88] (from the slanted part of the hysteresis loop) at R = 2.2 in system with N = 400 3 
spins, Ji n f = -k and J — 1. Again results agree with Eq. (j23H28|) with exponents given in 
Table □ Refer to Fig. Hi 

C. Finite Size Effects 

We study finite size effects in our simulation (mean field and IRM) by examining higher 
order spectra for various system sizes. We find that for smaller system sizes the high fre- 
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Vl.5 


v 2 


Ql.5 


Q2 


V c 


MFT 


2 


2 


3 


5 


2 


MF Sim. 


1.93 ±0.10 


1.92 ±0.12 


2.95 ±0.10 


4.93 ±0.12 


1.92 ±0.15 


IRM (d = 3) 


1.80 ± 0.07 


1.80 ±0.05 


2.73 ±0.06 


4.46 ±0.07 


1.82 ±0.14 


Experiment 


0.93 ± 0.20 


0.66 ±0.12 


2.81±0.08(h.f.) 
1.66 ± 0.12(l.f) 


4.39±0.15(h.f.) 
2.48±0.10(l.f.) 


0.63 ± 0.06 



TABLE I: We present the values of the exponents: V1.5, V2, Q1.5, and Q2 given in Eqs. (|23I27|) 
for mean field theory (MFT), mean field simulation (MF Sim.), the infinite range model (IRM) 
for d = 3, and experiment. While the exponents values for MFT were determined analytically, 
the exponents for the MF Sim., IRM, and experiment were determined through a non- linear curve 
fitting of the data. In particular, in our experimental plots (Figs. OHJ) we use the following 
window sizes to fit the exponents: fx : [1kHz, 200kHz] for Vx.5, fx : [20Hz,2kHz] for V2, fx ■ 
[14kHz, 40kHz] for Qx.5 (high frequency exponent), fx : [.4kHz, 18kHz] for Q2 (high frequency 
exponent), and fx ■ [20Hz, 2kHz] for Vc- We find that that these exponents do not change (outside 
of their error bars stated above) when the window size for their measurement is changed within the 
scaling regime of the data. Also, l.f. stands for low frequency and h.f. stands for high frequency, 
since for the experiment there are distinct l.f. and h.f. exponents, in some cases. 

quency scaling (and flattening due to the background term) is unchanged for second spectra, 
real 1.5 spectra, and real cross second spectra. However, at low frequency the scaling regime 
of the second spectra, real 1.5 spectra, and real cross second spectra rolls over (in the IRM 
and MF simulations) at frequency fx — jr—, where T max is the maximum avalanche dura- 
tion. T max is system size dependent where T max ~ L z , and L = N l l d . InaiV = 400 3 system 
(IRM) we find that T max ~ 4300. 

VI. DISCUSSION 

The mean field theoretical calculation was for a single infinite avalanche while the 
mean field simulation was obtained from a train of avalanches. Consequently, the train of 
avalanches introduces intermittency that lowers the magnitude of the mean field exponents, 
since the intermittency effectively adds white inter-avalanche noise to the intra-avalanche 
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noise seen in the MF calculation. From Table H] we notice that the exponents for our mean 
field simulation are systematically smaller (by an amount of 1% to 4%) than the expo- 
nents determined in mean field theory. Nevertheless, despite this small deviate, our mean 
field simulation results agree very well with mean field theory, corroborating our theoretical 
calculation. 

The background components (Q = 0), given by ri. 5 (/i) and T 2 (fb,r), are given in Fig. [T] 
for mean field theory and the mean field simulation. For the corresponding exponents Qi.5 
and Q 2 we find excellent agreement between mean field theory and the mean field simulation 
results. With the help of Uj we ascertain the following exponent relation for Q 2 and Q1.5: 

Q 2 = - 2 (29) 

avz 

Q1.5 = — + 1 (30) 
avz 

Plugging t and I /avz (given in Q]) into Eqs. (|29H30|) we find exact agreement in mean 
field theory: Q 2 = 5 and Q L5 = 3. For the IRM in d = 3 we find Q 2 = 4.40 ± 0.10 and 
Qi 5 = 2.72 ± 0.03, in close agreement with the table above. Also, for experiment we find 
reasonable agreement (within 10%) against the high frequency values for Q 2 and Q1.5; using 
Eqs. flUHHH and H we find Q 2 = 2.70 ± 0.05 and Q L5 = 4.02 ± 0.20 for experiment. Also, 
the measured experimental values for Q 2 and Q1.5 in the table above agree within error bars 
with the exponents for the IRM. 

Interestingly, we notice that the background components for experiment (given in Fig. |3J) 
have two scaling regimes: a flatter slope at low frequency, a steeper slope at high frequency, 
and a transition point at f cr oss — 320 Hz. This change in slope indicates that intra-avalanche 
correlations are effecting the power at f\ < fcross- 

The real 1.5 spectra and the second spectra are given in Figs. (|2l l3l6l l7|) . Since the 
background term is small for the real 1.5 spectra (for IRM, MF simulation, and experiment) 
we find that the real 1.5 spectra curves collapse upon themselves, in agreement with Eq. (|2*Bj) . 
The high frequency scaling exponent V1.5 for the real 1.5 spectra shows excellent agreement 
between mean field theory (V1.5 = 2) and mean field simulation (V1.5 = 1.93 ± 0.10). The 
second spectra (Fig. EJ) exhibits a conspicuous flattening due to the background term, except 
in the experimental second spectra (Fig. |7J) where the flattening in much less pronounced. In 
order to find the high frequency scaling exponent for the second spectra we do a non-linear 
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curve fit using the equation: AO * x~ A1 + (theoretical value, T ( 2 (/i)), where AO, and Al, 
are free parameters. We find very good agreement between mean field theory (V 2 = 2) and 
our mean field simulation (V 2 = 1.92 ± 0.12) for our second spectra exponent. 

The very weak dependence of the real 1.5 spectra on f\ and the strong fall off of the 
second spectra at high frequency suggests that the high frequency power comes from the 
fine structure of large avalanches (T >> 1/fi, T is the duration) and not small individual 
pulses (T ~ l//i)(5j. Through simulation we verify this claim in the mean field simulation 
and IRM. When we subtract all avalanches smaller than T = T max /4 from the Barkhausen 
train (where T max is the duration of the largest avalanche) and then determine the second 
spectra, we find no change in V%. However, when we subtract all avalanches larger than 
T = T max /4 from the Barkhausen train, we find that the second spectrum flattens and 
that there is an evident separation between the real 1.5 spectra curves (i.e. increased j\ 
dependence). In experiment we also find a weak dependence of the real 1.5 spectra on fi 
and a strong fall off of the second spectra, however, the second spectra for experiment falls 
off with an exponent V2 = 0.66 ± 0.12 versus an exponent of V2 = 1.80 ± 0.05 for IRM and 
Vi — 2 for mean field theory. This suggests that the high frequency power in experiment 
comes from the fine structure of larger pulses to a lesser degree than in IRM or mean field 
theory. This may be the result of dipole-dipole interactions that are present in experiment. 
Since the dipole-dipole interactions decay as a power law they may still be significant at 
short length scales, thereby resulting in suppressed spin flips that cause otherwise larger 
avalanches to be broken down in to smaller high frequency pulses. We are currently testing 
this hypothesis. 

In Fig. |U we give the real cross second spectra (r = f a /fb = 5) ^ or ^ e mean field 
simulation and the IRM, and in Fig. |S]we give the real cross second spectra for experiment 
(also r = |). The cross second spectra plots for r < | are similar to r = |, so to avoid 
redundancy r < ^ plots have been left out of the paper. The real cross second spectra not 
only strongly resemble the second spectra, but we also notice in Table H] that the exponents 
values for V 2 and Vc (where Vc was determined using the same non-linear curve fit used to 
find V 2 ) in mean field theory, the mean field simulation, the IRM, and in the experiment are 
identical or nearly identical. For the class of models and systems we study in this paper the 
real cross second spectra gives no new information, however, the real cross second spectra 
is useful when studying systems that have different dynamics on different length and time 
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scales fid]. 

Remarkably, we have found that Q1.5 and Q2 agree (within error bars) for the IRM and 
experiment. Since Q1.5 and Q2 are directly related to known exponents, as we have found 
above, they may be obtained from a standard analysis of power spectra (P(w) ~ w~~^) and 
the avalanche size distribution (D(S) ~ S~ T ) |2j. However, when we compare V1.5, V2, and 

between the IRM and experiment we find a significant difference, whose origin we are 
currently investigating. Thus by utilizing higher order spectra we present a more rigorous 
test for avalanche models against experiment. 

Our study of higher order spectra is a powerful tool to further our understanding of noise 
in disordered systems. In addition, our mean field results are applicable to a large array of 
systems, in particular systems discussed in Q]. 
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APPENDIX A: CALCULATION OF CORRELATION FUNCTIONS 
From Eq. (JIJ one can verify the following recursion relations: 



<n j > = < rij-i > (Al) 

< n] > = < n]_ x > + < rij-! > (A2) 

< n 3 j > = < n|_! > +3 < n|_! > + < nj-i > (A3) 

< n 4 . > = < > +6 < > +7 < n*_ x > + < %_i > (A4) 
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Using the fact that n a = 1 we determine the following initial conditions: 



< m > = 1 (A5) 

< nj > = 2 (A6) 

< n \ > = 5 (A7) 

< n * > = 15 (A8) 



Using the recursion relations and initial conditions we can find the explicit functionality 

of < n 2 >, < rij >, and < >: 



< n 2 > = < n 2 _ x > + < rij-i > 

— < n )-2 > + < %-2 > + < %-l > 

i-i 

... = < n\ > + 53 < n 3-m > 
m=l 

3-1 

= 2+Ei 

m=l 

= 3 + 1 

Thus we verify Eq. (jlj): 

< rtjUj+k >=< n 2 >= j + 1 (A9) 

Now let us consider: 



< ?^ > = < > +3 < n]_ x > + < > 

= < n]_ 2 > +3(< n 2 _ 2 > + < ri|_! >)+ < n^x > + < ra,-_ 2 > 

. . . = < n? > +3 ^ < n|_ m > + E < rij-m > 

m=l m=l 

= 5 + 3 53 (j - m + 1) + 53 1 

m=l m=l 

= ^(3j 2 + 5j + 2) 

Thus we find: 



14 



Now let us find < >: 



<n 3 >=i(3j 2 + 5j + 2) 



< n 4 j > = < n)_ x > +6 < r?j_ x > +7 < n 2 ^ > + < rij 

3-1 

= 15 + J2 I 6 < n )~m > +7 < n)_ m > + < 

m=l 

= I(6j 3 + 13j 2 + 9j + 2) 

So we have: 

<nJ>=i(6j 3 + 13j 2 + 9j + 2) 
Now let us look at the case where I > k > 0: 

< njn j+k n j+ i > = < rnn 2 +k > 

= < ^K-+fe-i + ^ 2 +fc _i) > 

= < n) > + < njn) +k _ x > 
= 2 < n 2 > + < njn 2 j+k _ 2 > 

= k <n 2 > + < n 3 > 

Thus: 

< njnj + krij + i >=< n 3 > +k < n 2 > 
Now let's look at the most general situation where m > I > k > 0: 



< njnj + knj+irij +m > — < njnj + knj +l > 

= < rijUj+ktri 2 ^^ + n j+ i_i) > 
= < n jn 2 +k > + < njn^ k n 2 Hl _ x > 
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< njrtj +k > + < n jnj+k (n 2 j+l _ 2 + nj-i- 2 ) > 
2 < n jn 2 j+k > + < njn j+k n 2 j+l _ 2 > 



(l-k) < n jn ) +k > + < n jn ] +k > 

(I - fc)[< n 3 > +k < n 2 >]+ < TijT& +k > 



Now let's determine < njn 3 +k >: 



< n jn 3 j+k > = < njin^-L + 3n) +k _ x + n j+k ^) > 

= < n i n| +fc _ 1 > +3 < rijU 2 ^^ > + < njn j+k -i > 

jfc-i 



< rij > +3 ^ < njn 2 j+m > +k < n 2 > 



m=l 
fc-1 



< > +3 ^ [< > +m < >] + k < n 2 > 



m=l 

k 

= < n 4 > +3k <n 3 > +-(3k - 1) < n 2 > 
j 2 

Thus combining the above results we verify Eq. (JJJ): 

k 
2 



< rijUj+kTij+inj+m >=< n 4 > +(2k + I) < n 3 > +-(k + 2/ - 1) < n 2 > (A13) 



APPENDIX B: TIME ORDERED PRODUCTS 

In order to evaluate Eq. we need to consider the time ordering of the indices of n, 
consider the following (i < N): 

y = n x + n 2 + ■ ■ ■ + njv-i + n N (Bl) 
Xi = Hi + n i+ i + . . . + n N -i + n N (B2) 

Now we want to evaluate y 2 , y 3 , y 4 we must write the expansion as a sum of time-ordered 
products, here's how we do it: 



r.2 



(ni + x 2 f = n\ + 2n x x 2 + x\ (B3) 
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xl = (n 2 + x 3 ) 2 = n 2 2 + 2n 2 x 3 + x 2 3 (B4) 
^3 = ("3 + Xif = n i + 2n 3 x 4 + x\ (B5) 



Thus by we can write Eq. (jB3|) as: 



(B6) 



iv 



f = EK 2 + 2^,,+ 1 ] (B7) 



i=l 

N N 

= £k 2 + 2n, £ %] (B8) 

1 = 1 J =2 + 1 

(B9) 

So we see that we have successfully written y 2 as a sum of time-ordered products. Now 
we can easily evaluate < y 2 > using Eqs. ©-(HI). This remarkable method can be used for 
higher powers to evaluate y 3 , and y 4 . 



APPENDIX C: CALCULATION OF HIGHER ORDER SPECTRA 



1. Second Spectra 

We first rewrite the second spectra as the follows: 



{FiiHfrNftFiiHfaN)}) 

D2{U!, 1\ ) — 



{H(i,N))i 
v-1-0 

S m 

(H(i,N))\ 



T/N-l T/N-l-6 

E e 2 ™ 6 " N ' T E (H(i, N)H(i + 8, N)) 

0=0 i=0 



Where % = t/N, (...)$ denotes sum over i (J2i=o ), an d Fi{..} represents the discrete 
Fourier transform over i. The variable u is conjugate variable to % in the Fourier transform; 
ji = uo/N is conjugate to t (the original time) since t = Ni. We will use to in subsequent 
calculations for consistency and do a change of variable to f% in our final expression. Also, 
we have set At = 1. 

Let's first consider the second spectra, when 9 > we have: 
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(H(i + e,N)H(i,N)) = ((xt +0 (N)-xr +0 (N)) 2 (x+(N)-xr(N)) 2 ) (CI) 
This product is already partially time ordered since we know 6 > 0, however when 8 = 0: 

(H(i,Nf) = ((xt(N)-xr( N ))±) (C2) 

In this case since 9 = no simplification can be made, there is no partial time ordering 
like in the case of Eq. (jCl|) . 

We begin by evaluating Eq. (jC2|) . which can be written as a sum of 4-pt functions 
(Eq.(jZJ)). To write Eq. ()C2|) as a sum of 4-pt functions we must use the method discussed in 
Appendix El which allows us to write any power of xf(N) as a sum of time ordered products 
in n: 



(H(i,Nf) 

= mf - 4(^) 3 sr + ^m^f - ^t(K? + (sr) 4 ) 



N/2-1 r 



N/2-1 



= E [4(iV/2 - 1 - j) + l] (nt N+j ) + 6 £ [2(JV/2-l-*) + l](4 +J f4 +k > 
j=o L k=j+i 

N/2-1 N/2-1 

+4 ]T [(3(N/2-l-k) + l)(n iN+j n 3 m+k )+3 E (2(N/2 -1 - I) + l)(n iN+ j n iN+kf^ N+ i) 

k=j+l 
N-l 

-2NY 



l=k+l 



N/2-1 



N/2-1 



(nf 



iN+j 



3=0 
N/2-1 N-1 



)+ E [ 3 ( n iN+j n m+i) + 3 ( n iN+jrf N +i) + 6 E ( n iN+jni N+ in 2 iN+m ) 



1=3+1 



m=l+l 



N/2-1 



+ 6 E E [ 2 ( N - 1 - A;) + l] [(n 2 iN+j n 2 iN+k ) + 2 ^ (n iN+j n iN+ in 2 iN+k ) 



j=0 k=N/2 l=j+l 
N/2-1 N-1 [ JV-1 

( 3 (iV - 1 - k) + l){n m+j n 3 iN+k ) + 3 ]T (2(JV - 1 - I) + l){n iN+j n iN+k n 2 N+l ) 



-4 E E 

j=0 k=N/2 I l=k+l 
N-1 I" N-1 

+ E [W - 1 - j) + l] (nf N+j ) + 6 E [2W2 - 1 - *) + l] (n 2 N+J n 2 N+k ) 

j=N/2 L k=j+l 
N-1 N-1 

+4 E [(3(^ - 1 - k) + l)(n iN+j n 3 iN+k ) + 3 E ( 2 (^ - 1 - I) + \){n lN+3 n lN+k n 2 N+l ) 

k=j+l l=k+l 

Performing the above sum we find: 
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(m,Nf) (C3) 
N [64 + 56(21 - 20i)iV 2 - 420iV 3 + 14(58 + 80z)iV 4 - 105iV 5 + 4(61 + 70i)iV 6 ] 



13440 

Finally summing over % we obtain: 



i(C4) 



T/N-l 

UN)= E (H^N) 2 ) (C5) 

i=0 

= ^N(N 2 + 2) 2 T 2 + 0{T) (C6) 

~ _AT 5 T 2 (C7) 
96 



Now let's consider the case where 6 > 0: 



{H{i + 0,N)H{i,N)) 

= - xr(JV)) 2 (x+^(JV) - x- +0 (N)Y) 

-N/2-1 N-l 

= E [2(JV/2 - + E [2(iV -1-0 + 1 

- (=0 l=N/2 
N/2-1 N/2-1 N/2-1 N-l 

E E ( n iN+jniN+kn 2 ( i+ g)N+l) + ( n 'iN+j n2 (i+9)N+l) E E { n iN+j n iN+kn 2 i+e ^ N+l ) 

j=0 k=j+l J j=0 fc=Ar/2 

JV-1 N-l 

+ E [ 2 E ( n iN+j , niN+kn' 2 i+e - )N+l } + (n 2 iN+j n 2 i+e)N+l ) 
j=N/2 k=j+i 

= —N 3 {N 2 + 2)(5iV 2 + 4 + 4i(iV 2 + 2)) 
576 

After summing over i and performing Fourier transform we find: 



T/N-l 

E • 



T/N-l-6 



E {Hit + 9, N)H(t, N)) 



t=o 



N 2 {N 2 + 2) 2 T 



From Eq.(H5J we know (H(i,N)) 2 = j^(N 2 + 2) 2 , also f 2 = u/N, and = 1/JV, so 
finally we obtain: 



S2V2J: 



1) - ™ 2 
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2. 1.5 Spectra 



We first rewrite the 1.5 spectra equation in a form we may analytically evaluate: 



Si. 5 (u,N) 



[Fj{v(i, N)}Fj{H(i, N)}) 
(H(i,N))i 



T/N-l T/N-l-e 

E e 2me " N / T E (v(i,N)H(i + 6,N)) 
e=o i=o 

{H(i,N))i 



Where v{i, N) = xf(N) + x ti (N). Now we consider the 1.5 spectra, when 9 = we have: 



T/N-l N/2-1 N-l 

ri. 5 (iV)= E (v(t,N)H(i,N))=[ E E "iiv+J 

i=0 £=0 i=AT/2 



(xt(N)-x-(N)y 



T/N-l 



= E 



i=0 



N/2-1 



N/2-1 



E [(3(7V/2-l-i) + l)(n^ +j )+ E mN/2-l-k) + l){n* N+j n iN+k ) 

3=0 k=j+l 



N/2-1 N/2-1 N/2-1 N-l 

-N/2 E [(nhf +j ) + 2 E (n lN+J n 2 lN+k )] - E E [ 2 (^ " 1 " k )( n iN +j 

3=0 k=j+l j=0 k=N/2 

= ^N{N 2 + 2)T 2 + 0{T) 

~ J- iV 3 T 2 
24 

For the case 9 > we find: 



71 



iiV+fc/ 



AT-1 



N)H(i + 9, N)) = E n iN+l (xt +e (N) - x~ +e (N)f 

1=0 



N-l 

= E 

1=0 
N-l 



N/2-1 



E [(2(N/2 -l-j) + l-N] (n iN+l n 2 i+e)N+j ) 



3=0 



+ E [2(N-l-3) + l\(n tN+l n% +e)N+J ) 

j=N/2 

= ^N{N 2 + 2){2Ni + N + 1) 
After summing over i, performing Fourier transform, and taking the real part we find: 



T/N-l 



Re{ E 



^2-kiBw/T 



T/N-l-d 



E (v(i,N)H(i + 9,N))} 



i=0 
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N 2 (N 2 + 2)T 



12uj 2 



Now we normalize the result with (H(i,N)) t = -^(N 2 + 2), and substitute f 2 = uj/N, 



and h = 1/N. 



ite{5i. 6 (/2,/i)} = ^ + ^T (C9) 



3. Cross Second Spectra 

We rewrite the cross second spectra as follows: 



(H(i,N)) f (H(j,M)) 

T/N-l T/N-l-e Ni/M 



3 



E e 2 ^ N l T E f E {H(j,M)H(t + 6,N)) 



0=0 i=0 N j=N(i-l)/M+l 



(H(i, N))i(H(j, M))j 

Where M = j^, N = -j-, and denotes sum over j (Ejicf" 1 )) Before we calculation 
of the cross second spectra we find it useful to define the following notation: 



M-l 

yj( M ) = J2 n Mj +i 

l=M/2 
M/2-1 

vti M ) = Yl n Mj+i 

1=0 

a = Ni/M 



b = N(2i + l)/(2M)-l 
c = N(i + 1)/M — 1 



c N-l 

x~(N)= £ (yUM)+y k (M))= £ n Ni+j 

k=b+l j=N/2 
b N/2-1 

xt(N) = J2(y+(M)+y k (M))= J2 n m+j 

k=a j=0 

Y k = yUM)+y^(M) 



X J =yl(M)-yj(M) 
Again first we consider case where 9 = 0, using the notation above we have: 
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T/N-l M c 
i-0 iV j=a 



T/N-l M c 

E Z E<Gf M - yj(M)) 2 (zt(N) - *r W) 2 > 



M 
iV 



=0 ^=a 

r/jv-i c 



E E*l[En- E n 



j=6+l 



T 2 (iV,r)= J] ^7 ^2(H(j, M)H(i, N)) (CIO) 

(Cll) 
(C12) 
(C13) 

Now we make the follow approximation for Eq. (jC13|) . we simplify the limits by setting 
N/M = 2 in the limits. By deduction this approximation will preserve the scaling for any 
N/M (multiples of 2) in our final answer. Using N/M = 2 to simplify the limits we find: 



i=0 j=a 

M T/N-l b c b c 



a = 2z 
6 = 2z 

c = 6 + 1 = 2i + 1 



We now expand Eq. ([013)1 : 

M T/N-l 

= ^ E (xl + Xl i+1 )(Y 2i -Y 2i+1 ) 2 

iV i=0 
M T I N - X 

\ ^ / y2 v 2 _L V 2 _L V 2 _L V 2 9 "!T 2 V V 9 "!T 2 V V \ 

~ JJ 2.^ \ Ji 2i Y 2i + A 2i+I r 2i+1 + Ji 2i+l Y 2i + A 2i r 2i+1 ~ ^2i Y 2i Y 2i+l ~ ^ 2i+l Y 2i+l Y 2i ) 

We evaluate the above terms below as a function of k, where k represents either 2i or 
2i + 1: 



Xl +1 Y k+1 Y k 

M-l M-l 

= E E 

m=0 q=M/2 I 
M-l 

+3 £ [2(M 1 ,5 1 ■ ■■ 2 

l=q+l 



[3(M - 1 - g) + l] (n k M+ m n 3 ik+1)M+q ) 
— 1 — + 1 ( n kM+ m n(k+i)M+ q nf k+1 - )M+ i) 
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M-l M-l M/2-1 



[2(M - 1 - /) + 1 (ntM+m1(fc+l)M+,«(fc + l)M+!/ 



-E E E 

m=0 i=M/2 9=0 

M-l M/2-1 M/2-1 

~ M E E [( n fcM+ m n^ +1)M+i ) +2 ^ (n fcM + m ^(fc+i)Af+/^(fe+i)M+ 9 ) 

m=0 «=0 g=2+l 
M-l M/2-1 r 

+ E E [ 3 ( M / 2 - 1 - q) + l] K« +m 4 + i)M +! ) 

m=0 q=0 L 

M/2-1 

+3 [ 2 ( M / 2 - 1 - + l] (nkM +m n { k+i)M+ q n 2 {k+l)M+l ) 

l=q+l 



Y 2 V 2 
uX k+l I k 

M/2-1 



M-l 



[2(M/2 - 1 - m) + 1 - M] + ]T [2(M - 1 - m) + 1 



rre=0 
M-l 



m=M/2 



X 



n=0 



M-l 

[ E (( n fcM+n n (fc+l)M+m) + 2 E ( n kM+nnkM+infk+l)M+n}) 

l=n+l 



X 2 k Y k 2 



[4(M-l-j) + l]K 



l Mk+ji 



M-l 

= E 

i=M/2 
M-l 

+6 ]T [2(M-l-m) + l](n Mfc+J n Mfc+m ; 



n=j+l 
M-l 



M-l 

+ E 4 [(3(M-l-m) + l)(n Mfc+J n Mfc+m ) 

rrj,=i 4-] 



l Mk+U 



m=j+l 

M-l 

+ E 3 ( 2 ( M - 1 - Z) + l) (n Mk+j n M k+ m n\ 

l=m+l 
M/2-1 M-l 

~ 2 E E [ 2 (M - 1 - m) + l] x [(n 2 2 

j=0 m=M/2 
M/2-1 r 

[4(M/2-l-j) + l](n 4 Mfc+J ) 



M/2-1 

Mk+j n Mk+m) + E 2 ( n Mfc+j n Mfc+Z^Mfc+r, 
i=J+l 



E 

j=o L" 

M/2-1 

+6 J! [2(M/2- 1 -m) + 1 

m=j+l 
M/2-l 



( n Mk+j n Mk+m) 

Mj 2—1 

+ ]T 4[(3(M/2-l-m) + l)(n Mfc+J 

m=j+l 



Mfc+ra/ 



23 



M/2-1 

+ Yl 3(2(M/2-l-/) + l)( 

n Mk+j n Mk+m n \lk+l) 

l=m+l 



kM+m) + 3 X! ( n kM+m n kM+l) 
l=m+l 

M-l 



X 2 k Y k Y k+1 

M-i r Af-i 
= M £ 

m=M/2 
M-l 

+3 X! [(^M+m^fcM+i) + 2 E ( n kM+m,nkM+inl M+q ) 
l=m+l q=l+l 
Af/2-1 M-l M-l 

E X! [("■fcAf+m^Lf+i) + 2 X! ( n kM+ m nkM+inl M+q ) 
m=0 l=M/2 q=l+l 

M/2-1 M-l M-l 

-m E E [P 

m=0 i=M/2 
M/2-1 

+M E 

m=0 L Z=m+1 
A//2-1 Af/2-1 

+3 E [( n fcM+m«Lf+«) + 2 E ( n kM+ m nkM+in 2 kM+q ) 



n kM+m n kM+l) + 2 E t \ n kM+mnkM+ q n 2 kM+ i) 
q=m+l 



(n 



M/2-1 

i \ i Q ^,2 



fcM+m 



)+3 E ( 
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n kM+l) 



kM+m lb kM+ 



l=m+l 



q=l+l 



X 



n kM+m n 



M-l 

E [ 2 (M-l-n) + l 

n=0 
M-l M-l 

E E 

n=0 m=M/2 
M/2-1 M-l 

~ 2 E E ) + 

m=0 i=M/2 

M-l 



M-l 



(fc+l)M+n) + 2 E ( n kM+mnkM+infk+l)M+r 
l=m+l 



M/2-1 
^— > 

m=0 

Now we can evaluate Eq. (jC13|) (where r = M/N): 



M/-Z-1 M-l 

E (( n fcM+m n (fe+l)M+n) + 2 E ( n fcA/+m^fcM+/^A: + l)M+n)) 
m=0 i=m+l 



4 



144 v- 1 m 3 < 
The 6 > case is readily evaluated as followed: 



25M 5 r 
144 



-T 
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(H(j, M)H(i + 9, N)) = ((# - y-) 2 (x+ +e - K+e?) 

M/2-1 M-1 

[2{N/2 -l-l) + l-N]+ [ 2 ( N ~ 1 " + 1 

L 1=0 l=M/2 
■M/2-1 M/2-1 

53 [ 2 ( n jM+ m njM+knf i+e ^ N+ i) + (n 2 M+m n 2 i+e ^ N+l ) 

. m=0 k=m+l 
M/2-1 M-1 

~ 2 X! X! ( n jM+ n njM+kn(i+8)N+i) 

m=0 k=M/2 
M-1 M-1 

+ X! ( n iM+ m njM+kn 2 i+ Q) N+l ) + (n 2 M+m n 2 i+d ^ N+l ) 

m=M/2 k-m+1 

After summing over z and performing Fourier transform and taking the real part we find: 



T/N-l 



Re{ £ 



2iri.eu)/T 



T/N-l-9 



J2 (H(j,M)H(t + e,N)}} 



9=1 i=0 

MN(M 2 + 2)(iV 2 + 2)T 
144cu 2 



FromEq.fUnD we know (H(i, N)) t (H(j, M)) t = ^(iV 2 +2)(M 2 +2), / 2 = w/JV, /„ = 1/N 
and fb = l/M finally we obtain: 



Re{S 2 (f 2 J a J b )} 



r 25r E 

Tf! + la~ 



(C15) 
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FIG. 1: We present T2(/i,r) (r = jj- = 1 corresponds to r2(/i), and we set fx = /& for r < 1) in 
the MF simulation and in the IRM. r2(/i,r) collapses for r < 1, in excellent agreement with Eq. 

For high frequency r 2 (/ x ,r) ~ /f Q2 , (see Tabled. Inset: r 1 . 5 (/ 1 ) for the MF and the IRM. 
For high frequency ri.s(/i) ~ fi 5 , see TableDJ 
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FIG. 2: We present the Re{Si.s(f2, fi)} in the MF simulation and in the IRM. At high frequency 
-Re{<Si.5(/2, /i)} ~ f2 Vl ' 5 (g° t° Tabled- There is not a visible flattening present due to r^(/i) 
(9 = 0) term, since the magnitude of this term is small relative to the / 2 dependent term. 
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FIG. 3: We present the S 2 (f2, h) in the MF and IRM. Notice the flattening due to the I2 (/i) 
{9 = 0) term. Inset: S^~(/2,/i) is the second spectra in the MF simulation and the IRM with 
background term subtracted. The bold lines adjacent to the MF curve is a power law with an 
exponent of -2, and the bold line adjacent to the IRM curve has an exponent of -1.8. 




FIG. 4: We present the Re{S 2 {f2, fa, fb)} for r = f a /f b = \ in the MF simulation and the IRM. 
Notice the flattening due to the r^^/^r = \) (9 = 0) term. 
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FIG. 5: We present T2(/i,r) (where r = = 1 corresponds to r2(/i), and we set /i = /& for 
r < 1) in experiment. r2(/i,r) collapses for r < 1, in excellent agreement with Eq. (|25|). For 
high frequency T2(/i,r) ~ fi® 2 , see Tabled Inset: ri.5(/i) for experiment. For high frequency 
ri.5(/i)~/f Ql - 5 , see Tabled 




FIG. 6: We present the real 1.5 spectra i?e{5i.s(/2, /i)} in experiment. At high frequency 
Re{Sl.5(f2, /l)} ~ /2 Vl ' 5 (g° to Tabled). The high frequency scaling regime is small for the 
experimental real 1.5 spectra since curve rolls over at f cross — 320 Hz. 
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FIG. 7: We present the second spectra, S , 2(/2,/i), in experiment. At high frequency the scaling 
is given S2(/2, /i) ~ f^ 2 , where V2 = 0.66 ± 0.12, that is significantly smaller than the mean 
field simulation (V 2 = 1.92 ± 0.12) and the IRM (V 2 = 1.80 ± 0.05). Further more, there is no 
conspicuous flattening in these experimental second spectra curves, indicating that the background 
term, T 2 N \fi), is small for experiment versus the mean field simulation and the IRM. Also, there 
is a noticeable separation between curves at low frequency that is not present in the mean field 
simulation or IRM second spectra curves. Inset: S^C.fej/i) is the second spectra in experiment 
with the background term subtracted, as a result the seperation between the curves vanishes, 
within error bars. 
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FIG. 8: We present the real cross second spectra, Re{S2{f2, fa, fb)}, for r = f a /fb = \, in experi- 
ment. Again, as in the experimental second spectra there is no conspicuous flattening due to the 
background term, T^\f^,r = ^). However, the seperation between the curves vanishes, within 
error bars, when the background term is subtrated off. 
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